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Abstract 

This work presents a new evolutionary optimization algorithm in the- 
oretical mathematics with important applications in scientific computing. 
The use of the evolutionary algorithm is justified by the difficulty of the 
study of the parameterization of an algebraic variety, an important prob- 
lem in algebraic geometry. We illustrate an application, Evo-Runge- 
KuTTA, in a problem of numerical analysis. Results show the design 
and the optimization of particular algebraic variety, the explicit s levels 
Runge-Kutta methods of order q. The mapping between algebraic geom- 
etry and evolutionary optimization is direct, and we expect that many 
open problems will be modelled in the same way. 

1 Introduction 

In science and engineering, problems involving the solution of a polynomial 
system happear everywhere. In this work we study the optimization of a positive 
real values function Q defined over X C M^, the set of real solutions of a system 
of m polynomial equations /i(x) = for i — I, ... ,171. 

The key problem of optimizing G over X is the necessity to control X via 
its parametrizzations: for reasonably large values of N and m, it is an open 
problem to find the connected components {Ui} of X , their local dimensions di 
and their local parameterizations {di : R''* Ui}. Evidences does not suggest 
the use of symbolic computation or numerical method [24] . For this reason, we 
state a particular methodology based on Evolutionary Algorithms. 

The Evolutionary Algorithms (EAs) are a generic population-based meta- 
heuristic optimization algorithm bioinspired [18]. They evolve a population of 
candidate solutions for the problem by reproducing and selecting with respect 
to a fitness function. They recently improve on several field because it is not ne- 
cessary any assumption about the solutions space of the problem. In particular 
the evolutionary algorithms fit very well into the optimization problem where 
the fitness function is cheap. In this research field, the EAs were used to solve a 
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particular class of polynomial system in [26] and an hybrid algorithm was pre- 
sented with the same aim in [28]. In addition there were genetical approaches 
in pure algebra in [24] and [4]. 

In applications it is not usual have a proper justification for the use of such 
optimization algorithms, apart cases where the dimension of the problem is 
prohibitive for any numerical computation. We show that algebraic geometry 
provides a justification for why it is important to use evolutionary optimization 
algorithm to optimize any positive real function over the algebraic variety X, 
but this justification holds for any kind of optimization over X, [3] [8], [12]. 

In Section 2 we introduce some notions about the algebraic varieties for a 
non expert reader. In Section 3, we set up all the theoretical objects used in 
the most general enviroment and after that we present the abstract algorithm 
Evo-Alg- Variety. 

We apply our suggested methodology in a numerical analysis interesting 
problem: in [1] , [22] , one wants to perform some particular features of the RK 
methods. With Evo-Runge-Kutta, a new evolutionary optimization algo- 
rithm designing Runge Kutta methods, we find new Runge-Kutta methods of 
order q with minimal approximation error. 

At first, we introduce the Rungc-Kutta methods in Section 4, sketching some 
important results in [19]. In Section 5, we apply them to define the problem 
from a geometrical point of view. We show that 

= {s-level Runge-Kutta methods of order accuracy g} 

is an algebraic variety and one cannot compute its dimension and its local 
parameterizations. Section 7 shows the obtained results. 

The final section discusses the conclusions of the work underlining the pos- 
sible directions of future optimization algorithms and the implications for alge- 
braic varieties and pure algebra. 

2 Using Algebraic Varieties 

In this section we sketch some notions about the algebraic varieties for the reader 
without this background; we suggesting to an expert in algebraic geometry to 
skip this section. 

Claiming a zero set of a polynomial system being an algebraic varieties has 
some underlying effects. One of the main differences conccirns the topology. 
When we write M*^, we normally mean the set of A;-tuples of real numbers with 
the standard (or Euclidean) topology, i.e. the topology in which we consider 
the open sets as the the fc-dimcnsional open ball, B{xq, e) = {x G R'^ : |x — 
a;o| < e}, of center xq € K*^ and radius e > 0. Instead, A'°(R) is the affine 
fc-dimensional space: the set of M-points is the same as before, but the topology 
is radically different; here, as in any algebraic variety, it is used the Zariski 
topology: a closed subset (the complement of open subset) is the zero set of 
some polynomials in k variables. In contrast to the standard topology, the 
Zariski topology is not Hausdorff^ . This non-trivial example shows the difference 
between the two topologies. 

topological space Y is Hausdorff if it is possible to separate two distinct points with 
two disjoint open sets. 
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Since the objects (the varieties) are different, so the maps among those wih 
be. Let us start, as previous, with the standard case and let f{x) = i be a 
real function on an open subset of M; the function is not defined for x = 
and, similarly, the zero value is not in the image of /. This happens because 
the set {0} is a closed point. In the Zariski topology there are a lot of closed 
sets more; in particular algebraic map can be undefined over a finite number of 
those. For the reader who meets algebraic maps for the first time, it is useful to 
know that such maps have not uniquely way to be written. For instance, let X 
be the algebraic variety given by = + x"^ in A^(M) and let / : X — >^ M be 
f{x,y) — {y/x)"^] it seems that / is not defined in a; = 0, but using the defining 
polynomial of X, one has 

f{x,y) = [-) = — = 2 — =x + l. 

\XJ X'' X'' 

So, for X = 0, f assumes value 1 and / can be written either as {v/x)^ or as 
a; + 1. I refer to [21] for the undefined geometrical notions in this section and 
in the nexts sections. 

2.1 Local Dimension and Local Parametrization 

The first natural question regarding the geometry of an algebraic variety V is 
about its dimension, or, with different words, the number of free parameters 
that we need to use to parametrize it, globally or locally. 

A very nice introduction to dimension theory is in Chapter II. 8 of [13], in 
which there are equivalent definitions for the dimension of an algebraic variety: 
for example using the Noether Normalization Lemma (sec [2]). However, there 
exists a finite^map (p : V ^ A'^(M), where d is the dimension of V. It is 
extremely difficult to compute d (see [7], [11]): there are some methods in 
computational algebra where the complexity of the computation depends on 
the number of generators m, the number of variables N and the degree of the 
polynomials q as m°^^^ q°'^'^^ (see [16], [17]). Hence, it is clear that with A'' and 
m bigger enough any symbolic approach to this problem fails. 

Moreover, since we work with real numbers, there is an additional compli- 
cation: working with the complex field guarantees that the varieties are con- 
nected. For instance the circle Cc given by -I- = — 1 is a one dimensional 
variety in the affine complex plane and, if we compactify it, one knows that 
Cc is projectively equivalent to (the compactification of) the hyperbola 2c, de- 
fined by xy = 1. Instead, working in the real field, Cr is the empty set and 
2r is non-connected: Ir is a union of two pieces Ui = {xy — l,x,y > 0} and 
U2 = {xy = l,a;,y < 0}. Trivially, Cr and Ir are far to be equivalent, just as 
sets. 

Thus, an affine real variety V should be seen as a union of connected com- 
ponents, V — Di^iUi] each one of this component Ui has a proper dimension di 
with 0<di<d. 

^An open subset (7 of a scheme X is affine if there exists a ring R such that U = Spec{R); 
the Zariski topology, that we have discussed before, comes from the notion of Spec of a ring. 
A map between algebraic schemes, f : X Y, is 3. finite map if there exists an open aiRne 
covering of Y, {Ui = Spec(Bj)}ig/, such that for each i, f~^{Ui) = Spec{Ai) where Aj is a 
finite Bj-algebra. 
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In order to perform any kind of optimization over an algebraic variety V, 
we need to parametrize it. In fact, we need to control the set V using free 
parameters {pi ,.. .,pk} where k is at least the dimension d of V. More explicitly, 
given a connected component decomposition, V = Ui^jUi, we need to know the 
functions 

for all i £ I. This problem, translated in the algebraic geometry form, 

is the problem of finding the local parameterization of an real algebraic variety. 
1 remark, that, even if the meaning and the appearance of the two 9i are exactly 
the same, the algebraic one carries all the differences discussed. 

If we narrow down our investigation, for a moment, and consider only curves 
(one dimension algebraic varieties) in the afRne plane (see [15], [25]), the que- 
stion is the following: is the curve rational? We know that non-rational curves 
exist. So not all varieties X are birational equivalent (generalization of rational 
concept) to P''(]R). The theory of Groebner basis (see [9]) or the Newton Polygon 
(sec for recently connection with tropical geometry [10]) could be applied to 
tackle this problem, but, with a lot of generators, the computational time is not 
acceptable. 

The problem is complicated, hence, excluding some particular case, finding 
local parameterizations of an algebraic variety generated by a large number of 
polynomials in a large afHne space is an open problem (see [29]). 

3 The General Algorithm for Optimization over 
an Algebraic Variety 

In this section, we state the optimization problem and the suggested algorithm 
in the most abstract description. 

Let X C be the set of real solutions of a system of m polynomial equa- 
tions 

' /i=0, 

/2=0, 

< . 

. fm = 0. 

Let N and m be great enough. Let Q he a, positive real values function defined 
over X and let us consider the optimization problem consisting of finding x € X 
minimizing the value ^(x) in G{X) C M+. 

In applications (see for examples [23], [27]), it is usual subdividing the gen- 
erators in a finite number of sets. Let {fi}i=kj-i,...,kj be this subdivision in c 
sets (with fco = 1 and kc = m). We define the algebraic varieties 

X, = {xeA^ : /.(x)=0, V* = l, ...,%}, 

for j = 1, . . . , c and Xq = Aj^ and Xc — X. Hence, there is the following chain 
of c inclusions: 

A]^ = ^0 3 Xi 3 • • • D Xc-i 3 Xc = X; 
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if {/i}i=i,...,m is a minimal basis for the ideal (see [9]) of X, then inclusions are 
strict. 

If there is no motivation to make a subdivision of the generators, then it is 

possible to choose, ior < i < v, ki ~ i[™/t)J, with v a suitable positive integer: 
the length of the chain will be exactly c = v + 1. One should choose v coherently 
with the hardness of solving the system of polynomials. 

Applications (sec [5]) also suggest that some constrains fi{x) = arc more 
important than the others. So, let r = (ri)i=i,...,m be an element of the simplex 

m 

Am = {(ri)i=i,...,„ : n > 0,^ri = 1}; 

i=l 

Vi is the weight of the polynomial equation fi{x) = 0. Thus, for j = 1, . . . , c 
the positive real value function Hj : — >■ M_|_ is defines as 

^.•(^) = Eni/i(x)|. 

i=l 

Let e = {ei < e2 < ■ • ■ < Cc} be a set of positive non increasing real numbers; 
we define the e-tubular neighbourhood of Xj being 

e{Xj) = {x€A^ : Hk{x) < efe, fc = 1, . . . , j}. 

Finally, we propose the algorithm in Table 1 and we describe the objects 
used in Table 2. 

Evo-ALG-VARiETY(iV, {/,}, {%}, Q, r, e, {Sj}) 

1. 5*0 = 0; 

2. For {y from to c — 1) 

3. EA-minimize {T-Ly+i over e{Xy) starting from Sy); 

4. save-good-solution (Sy+i); 

5. End for; 

6. EA-minimize {G over e{X) starting from 5c); 

Table 1: The Algorithm for the optimization of Q over the algebraic variety X 

The kernel of the optimization is the fact that we sequentially produce an 
homogeneous covering of the affine varieties Xj: in each y-cycle we optimize 
Hj+i on e{Xj). In other words, the fitness function of every y-cycle coincides 
with the constrain of e{Xj+i). 

For this reason, in addition, we save the points x such that Hj+i {x) < ej+i , 
i.e we save the elements x G ^{Xj) such that x £ e{Xj+i). Those feasible 
solutions are the suggested starting points for the evolutionary algorithm of the 
next cycle. 

These information can be used for any optimization over the variety X . In 
fact, in the end of the algorithm, we run the last optimizzation over the feasible 
points of e{Xc) = e{X) minimizing the function G- 

Since we work in a real field and since the N and m are large enough to make 
impossible the analytic computation of the local parameterizations 6i : M.'^' Ui 
(with X = Ui^iUi), then the use of the EAs is fully justified. 

What fallows it is devoted to an application of this abstract algorithm to a 
numerical analysis open problem. 
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Object 



Description 



N 



number of variables 

number of equations 

polynomial equations 

index of the subdivision 

number of the sets subdivision 

function to optimize 

weights for the polynomial equations 

numbers definining the e- tubular neighborhood e{Xj) of Xj 
set of feasible solutions for e{Xj) 



m 



{/J 



c 



Q 



r 



e 



Table 2: Objects used in the algorithm desctibed in Table 1 



4 Runge-Kutta Methods 



Problems involving ordinary differential equations (ODEs) can be always refor- 
mulated as set of TV coupled first-order differential equations for the functions 
yi, having the general form 



A problem involving ODEs is completely specified by its equations and by 
boundary conditions of the given problem. Boundary conditions are algebraic 
conditions: they divide into two classes, initial value problems and two-point 
boundary value problems. In this work will consider the initial value problems, 
where all the are given at some starting value to, and it is desired to find the 
2/i's at some final point tf. In general, it is the nature of the boundary condi- 
tions that determines which numerical methods to use. For instance the basic 
idea of the Euler's method is to rewrite the dy's and dt's in (1) as finite steps 
6y and 6t, and multiply the equations by 6t. This produces algebraic formulas 
Sy with the independent variable t increased by one stepsize St; for very small 
stepsize a good approximation of the differential equation is achieved. 

The Runge-Kutta method is a practical numerical method for solving initial 
value problems for ODEs [19]. Runge-Kutta methods propagate a numerical 
solution over an iV + 1-dimensional interval by combining the information from 
several Euler-style steps (each involving one evaluation of the right-hand g's), 
and then using the information obtained to match a Taylor series expansion 
up to some higher order. Runge-Kutta is usually the fastest method when 
evaluating gi is cheap and the accuracy requirement is not ultra-stringent. 

In what follows, we summarize some basic notions about Runge Kutta me- 
thods (RKms) and we remark the results that we are going to use. 

Without lost of generality we can consider, instead of (1), the autonomous 
system 



given by the function g' : O — >■ M" (with O being an open subset of M") such 
that the Cauchy problem makes sense. Moreover, we suppose that g satisfies 



dyijt) 
dt 



9i{t,yi, - ■ ■ ,yN) with i = 1,2, ... , 



N. 



(1) 




(2) 
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the hypothesis of Schwartz' Theorem for mixed derivatives of each order that 
we will use. 

The Runge-Kutta are a class of methods to approximate the exact solution 
of (2). The structure of s-level implicit RKm is the i-stage value 

s 

ki = g{yo + h^aijkj), 
for i = 1, . . . , s. The numerical solution yi e of the problem (2) is 

s 

yi = yo + h^Wiki] 

where h is the step size. In the implicit case it is necessary to invert (almost 
always numerically) the function g to reach the ki values. To avoid that, the 
explicit RKms are often used: aij = 0, Vi > j and so the fc, are obtained 
sequentially. Despite this nice feature, there are some problems (the stiffness 
problems) where only implicit RKms are required. 

The Butcher Tableau [6], in Tables 3, shows all the used parameters: the 
tiii and the a;.j arc real numbers and characterize a given method with respect 
another one. A RKm for a non-autonomous system involves also the q para- 
meters regulating the step size in the time axis. In the autonomous restriction 
such a parameters can be avoided because they fulfill the following autonomy 
conditions: 

i-l 

Cj = ^ aij, Vi = 1, . . . , s. (3) 

3=0 



C\ 


ai,i 


ai,2 


ai,3 


Ol.s-l 


ai.s 


C2 


a2,i 


a2,2 


02,3 


0'2,s-l 


a2,s 


C3 


(13,1 


013,2 


013,3 


13,8-1 


a3,s 






as-1,2 


O's-1,3 




^3-1,8 








as,3 


0's,s-l 








W2 


Ws 


Ws-1 


Ws 



Table 3: A Butcher Tableau for an s-level implicit Runge-Kutta method. 

The order of approximation of a RKm is the order of accuracy with respect 
to h of the local truncation error, 

0-1 (^) = y{to + h) -yi{h), 

measuring how much the numerical solution yi (h) matches with the exact solu- 
tion y{to + h). 

To understand the relation between its parameters and its order of approxi- 
mation, in [19] it is developed a combinatorial interpretation of the equations 

that the parameters of a RKm must satisfy (so called order condition equations): 
there exixts a bijection between the set of p-elementary differentials and the set. 
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Tp, of rooted labelled trees of order p. Moreover, the following theorem shows 
how to obtain the condition equations of order 5 of an s levels RKm. 

Theorem 4.1. A Runge-Kutta method defined by a Butcher Tableau is of order 
q if and only if 

X:«;,$,(t) = ^ (4) 

for all the rooted labelled trees t of order at most q. 

Using this result, in [14] the system of equations is computed symbolically. 
One knows how to obtain Tp and so the number of order condition equations. 
It blows up like the factorial with respect to the order q (see Table 4) : this fact 
plays a key role in the choice of the evolutionary algorithm for the solution of 
the problem. 



order 



equations 



123456 78 9 10 



1 1 2 4 9 20 48 115 286 719 



Table 4: Number of the order condition equations up to order 10. 



We did not write the explicit form of ^j(t) and j{t) beacuse these involve 
more notions about the rooted labelled trees (we refer to [19] for more dectails). 
What we need to know, here, is that $j(t) is written as a sum of aij's products 
and j(t) is an integer number caming from the tree's structure. 

For the reader do not use to this numerical analysis branch, we sketch some 
examples. 

Example. The parameters of a two levels explicit RKm have to satisfy 

wi + W2 = 1; 

1 

W2a2,i = -. 

For instance the Euler method has wi = 0, W2 = 1 and 02,1 = 1. 
Example. For an implicit three levels RKm of order three, the parameters fulfill 

3 



Yl ^j'^J'fc = 2' 
j,k=i 

3 1 

,fe,(=l 
3 1 

Wjaj^kak,i 



3,k,l=l 
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The next theorem clarifies how the Taylor series of the numerical and the 
exact solution are. In fact, the /i^-terms of their series are expressed using 
G'^{t){yo), the J-component of the elementary differential of g corresponding to 
the tree t evaluated at the point j/o- For every G'^t){yo) an integer coefficient 
a{t) appears naturally to weight its contribute in the Taylor sequence. 

Theorem 4.2. // the Runge-Kutta method is of order p and if g is (p + 1)- 
times continuously differentiable, we have 

y'(yo + f')-y( = j^^^ E «(t)e(t)G-^(i)(yo)+^(/*^+') (s) 

where 

s 

e{t) = l-^{t)J2^j<^jit) (6) 

is called the error coefficient of the tree t. 

The error coefficient e{t) expresses the difference of the Taylor series of the 
numerical and the exact solution at the term of the elementary differential given 
hyt. 

The optimization problem that we want to solve is minimizing the local 
approximation error of order q + 1 in the set of RKms of order q: in other 
words, after finding feasible RKms of order q, our aim is to optimize such errors 
e{t) for all the trees of order q+1. 



5 Runge-Kutta methods as an Algebraic Vari- 
ety 

In this section we present the problem of minimize the local error of a RKm as 
an optimization problem over an algebraic variety. We define 

Vg = {s-level Runge-Kutta methods with accuracy order q} 

and we call EV^ its subset having only the explicit ones. 

Since the Table 3, the RKms have + 2s free coefficients. The autonomous 
conditions (3) show that the parameters Ci are dependent by {ai,j}. Moreover, 
the order conditions in Theorem 4.1 are polynomials, hence is an affine 
algebraic variety in the affine real space A*^*+^^(M), minimally defined by the 
following polynomials in s(s + 1) variables: 

E ^^-^i W = G Ti U T2 U • • • U T,. 

Thus, we rewrite the Theorem 4.1 as: 

Theorem 5.1. Let x = {{aij)i<j^i<:s, {wi,W2, ■ ■ ■ ,Ws)), then 

X is the parameters set of an s levels RKm of order g <;4> x e V^" . 
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Proof. Considering x as the parameters of an s levels RKm of order q, then x 
satisfies the order condition equations (4) and then x belongs to Vg as a point 
in the afifine space A*(*+^) (R). Viceversa if the point x of A*(*+^^ (R) belongs to 
Vg, then X satisfies the equations (4), that are also the order condition equations 
of the s-level RKms of order q. □ 

Similarly the algebraic variety EV^ in A * 2^ ' (M), minimally defined by the 
same polynomial equations and by {aij = Vi > j}, is the variety of the explicit 
s levels RKms of order q; EVg is also a subvariety of Vg. The Corollary 5.1 
holds too. 

Corollary 5.1. Let y = ((ai,j)i<i<j<s, (w'i,W2, • • ■,Ws)), then 

y is the parameters set of an s levels explicit RKm of order q '^y G EVg. 

Avoiding repetitions of EV^ and Vg , every results that we state for EVg is 
true also for Vg. 

We know that the EVg is composite by connected component: EVg = 
Lii^iUi. One knows that for large value of s and q the computational time for 
computing the local parameterizations, 

is not acceptable. In fact a symbolic approach of the problem is not feasible 
(look at [14]). 

Since minimizing the local error over the variety EVg involves those param- 
eterizations, then the use of the EAs is legitimate. 



6 EVO-RUNGE-KUTTA 

The aim of this research work is to obtain new explicit or implicit Runge-Kutta 
methods of maximal order q that minimizes local errors (6) of order q + 1. 
Now, We state the optimization problem for EVg , because we will show the 
numerical results about EV^, EV^ and EVf, but everythings hold with the 
opportune modifications for Vg too. 
We denote, as in Corollary 5.1, 

s(s + l) 

X = ((aij)i<i<j<s, (wi, W2, . . . , Ws)) G R 2 , 

where {aij} and {wi} are the coefficients of the Butcher Tableau of an explicit 
Runge-Kutta method, x is a RKm of accuracy order q (i.e. it lies in EVg) if 
and only if it respects the following constrains: 

1 

J2 u^j'^jit) = -7^, Vt e Ti U T2 U • • • U T,. 
Moreover, we want that x minimizes the local errors: 

s 

e(i) = 1 - 7(t) ^j^ji*)' ^* e T,+i 
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The impossibility of analytically computing the functions 0^ ^ and the hard- 
ness of analyze it numerically gives us the motivation to use evolutionary algo- 
rithms to face this difficult global optimization problem (see [18]). Moreover, 
we produce a remarkable set of feasible solutions that could be used for different 
optimization problems over the algebraic varieties EVg and Vg-. for examples, 
minimizing the computational time or maximizing the convergence area Sa for 
the implicit RKms. 

Fixed the number of levels s and the order of accuracy q the feasible condi- 
tions for a RKm being of order p are given by 



a{t)\e{t)\ 



< Cp, Vp= 1,2,...,? 



(7) 



where the Cp values is showed in the Table 5. We require an error large at most 
1 • 10~^^ for each constrain, but we use the amplification factor 4 to widen the 

tubular neighborhood of the varieties EVg and to develop a greater number of 
acceptable solutions, keeping the diversity of the population. 



order 


Number of condition equations 


Cp 


1 


1 


1 *4* 10-1^ 


2 


2 


2*4* 10-1^ 


3 


4 


4*4* 10-1^ 


4 


8 


8*4* 10-1^ 


5 


17 


17*4*10-15 


6 


37 


37*4* 10-15 


7 


85 


85*4* 10-15 


8 


200 


200*4* 10-15 


9 


486 


486*4* 10-15 


10 


1205 


1205 *4* 10-15 



Table 5: Values of Cp for each order up to 10 



The fitness function is 



Summarizing, we narrow down the general algorithm Evo-Alg- Variety 
(see Section 3) fixing X = EVg and N = «(s+i)/2. The minimal basis of gene- 
rators for EVg is given by (6) and we subdivide them in q sets, following the 
subdivision in labelled rooted trees. We also use the combinatorial definition of 
the Runge-Kutta order conditions to set the weight: in fact 



is the weight of the equation l~j{t) Ej=i = with t G Tp, as suggested 

by (5). The final setup is the definition of the Cp = Cp, concordly with Table 5. 
In conclusion, "Hq = Tq and Q = J^q+i as in (8). 
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Evo-Runge-Kutta(s, Qstart, A, fi, {cp}p=i g, acc — mean, acc — best, 

Imaxi best, oldSolutions, others par's) 

1. output = true; 

2. q = Qstart; 

3. While (output) 

4. output=Evo-Runge-Kutta-Cycle (s, q,...); 

5. q=q+l; 

6. End While; 



Table 6: The Algorithm Evo-Runge-Kutta 

This particular choice of Hq and Q simplifies the the algorithm: Evo- 
RuNGE-KuTTA has no final optimization of Q, because G = T-Lq+i and the 
minimizing problem of Q is seen as another repetition of the main cycle. 

However, in the Runge Kutta case, we do not know the dimension of the 
varieties EV^. This approach allows as to solve this complication. In fact, for 
any choice of the number of levels s, we always run the Evo-Runge-Kutta 
with q — 2. The algorithm explores the space of solutions, Aj^, finding points 
in e{EV2); after that, under the constrain being in e{EV2), it finds feasible 
solutions in e{EV^); and so on until the variety e{EVgj^i) is empty. Evo- 
RUNGE-KUTTA cannot prove that e{EVq_^_i) is empty. What happens is that 
the algorithm minimizes automatically the local error of order g + 1 in e{EVq ) 
(i.e. the fitness function Q) and it fails to find new feasible solutions in e{EVgj^i). 

In Table 6 we show the structure of the algorithm: the main cycle is explained 
in Table 7. Table 8 shows the parameters of the evolutionary algorithm. 

We use the package CMA-ES (see [20]) for the evolutionary optimization. 

7 The Numerical Results for EVi, EV^ and EVl 

In this section we show the numerical results obtained with Evo-Runge-Kutta 
for the explicit case in the levels s = 3, 4, 5. Since [6], the variety EV^ is empty 
and so the optimization for s = 5 runs over the same equations of s = 4 but 
with more variables; hence EV^ should be dimensionally bigger than EV^. The 
Table 9 shows the feasible points founded during the runnings: it respects what 
expected. 

In addiction, we present in Table 19, and respectively in Table 20, some of 
the Runge Kutta elements in Pareto front of the final optimization in EV^ , and 
respectively EV^: the condition eo,o = e(io,o) is zero for both these computa- 
tions (we enumerate the trees and their local error as in [19]). 

To have a proper view of the results, we list the Runge Kutta methods 
suggested in the Table 20 with all the digits from the Table 10 to 18. 

8 Conclusion 

The designed and implemented evolutionary algorithm, Evo-Runge-Kutta, 
optimizes implicit or explicit Runge-Kutta methods in order to find the maximal 



12 



Evo-Runge-Kutta-Cycle(s, g, A, /i, {cp}p=i,..._q, acc — mean, acc — best, 
Imax, best, oldSolutions, others par's) 

1. t=0; 

2. inizializePopulation(pop^*^); /* random generation of RKs */ 

3. m\i\dlvz.Q{new Solution)-, /* new array of feasible solutions */ 

4. evaluationPopulation(pop^*^); /* evaluation of RK systems */ 

5. While ((t< /maa:)&&bestError< acc — best) 



6 . copy (pop^ , popMut^ ^ ) ; 

7. mutationOperator(popMut^*'', fj); 

8. isFcasiblc(popAf7/i^,*'', oldSolutions); 

9. evaluationPopulation(popMut^*''); 

10. pop^*'''^'=selection(pop^*\ popMut^^^); 

11. computeStatistics(pop^*'''^^); 

' 12. saveSolutions(newS'oZ Mtion); 

13. t=t+l; 

14. End While; 

16. If {newSolution is empty) Then 

17. Return false; 

18. Else If 

19. Return true; 

20. End If; 



Table 7: The main cycle in the algorithm Evo-RUNGE-KUTTA 



Parameter 


description 


Value 


Qstart 


starting value for the order the RK-methods 


fixed 


S 


levels of the RK-methods 


fixed 


A 


number of elements of pop. 


1000 




number of crossing over for pop. 


500 


Ijnax 


maximum iterations of the cycle 


10^ 


acc — mean 


acceptable value for fitness's mean of the pop. 


Cg+l 


acc — best 


acceptable value for fitness of the best element 


c<i+i/4 


output 


output of EVO-RUNGE-KUTTA-CYCLE 


true/false 


oldSolutions 


points in e{Xq) 




newSolutions 


points in e{Xq+i) 





Table 8: Parameters of used in evolutionary algorithm Evo-Runge-Kutta 



Order/ Level 


3 


4 


5 


2 


46 


108 


34 


3 


146 


197 


140 


4 





364 


932 


5 












Table 9: Feasible solutions for the varieties EVg with s = 2, 3, 4 and g = 2, 3, 4, 5 
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0.6284799329301066 

0.0768995143272878 0.2946205527425946 

0.9213018835278037 -0.3307380161751250 0.4498926859290699 

0.1599508127125725 0.2455227275515280 0.4335530269375665 0.1609734327983330 

Table 10: The Butcher tableau of first 4 levels Rungc-Kutta in Table 20 

0.7363018963226102 

-0.0824111295277597 0.3461092332051678 

0.8341073930886564 -0.2376248889779167 0.4318380713817837 

0.1571049761428195 0.2866708210526255 0.3980679473446672 0.1581562554598879 

Table 11: The Butcher tableau of second 4 levels Runge-Kutta in Table 20 

0.7066207865932520 

-0.1986791412551434 0.4920583546618451 

0.9521842711595797 -0.3227487964624350 0.4140057189110243 

0.1495805086972162 0.2852315575222697 0.4127591061897394 0.1524288275907748 

Table 12: The Butcher tableau of third 4 levels Runge-Kutta in Table 20 

0.7559397531895379 

0.0103807700662207 0.2336794767442519 

0.7531440557235027 -0.1855193036061356 0.4515825081294965 

0.1620639462265475 0.2833431842164397 0.3921840325526135 0.1624088370043993 

Table 13: The Butcher tableau of 4-th 4 levels Runge-Kutta in Table 20 

0.7963012365770038 

-0.3628863293699752 0.5665850927929587 

0.9145513895258300 -0.2798774916728459 0.3940759984502829 

0.1451855083814838 0.3164608254865673 0.3906989174541415 0.1476547486778075 

Table 14: The Butcher tableau of 5-th 4 levels Runge-Kutta in Table 20 

0.7573881838384795 

-0.1741476827373300 0.4167594988988161 

0.8573072494974391 -0.2469806098316732 0.4182913673468296 

0.1534104115563865 0.2981759720267748 0.3935172041716111 0.1548964122452277 

Table 15: The Butcher tableau of 6-th 4 levels Runge-Kutta in Table 20 

0.6144690308782115 

0.0769053801528305 0.3086255889689522 

0.9475038121348336 -0.3525125955836914 0.4489168429216662 

0.1594196652643455 0.2418674636726206 0.4381018921031239 0.1606109789599099 

Table 16: The Butcher tableau of 7-th 4 levels Runge-Kutta in Table 20 

0.6957301007596283 

0.0048639630817450 0.2994059361585764 

0.8469017267939670 -0.2590107629375511 0.4434995354305374 

0.1594489916633499 0.2691847382771309 0.4110483067778322 0.1603179632816870 

Table 17: The Butcher tableau of 8-th 4 levels Runge-Kutta in Table 20 
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0.6313830140634179 

-0.0266776618862919 0.3952946478228191 

0.9776608341341720 -0.3627250371329542 0.4349091893706292 

0.1552061988645526 0.2544888797173599 0.4329627509848846 0.1573421704332030 

Table 18: The Butcher tableau of 9-th 4 levels Runge-Kutta in Table 20 



order of accuracy and to minimize theirs local error in the next order. The 
results presented in this article suggest that further work in this research field 
will advance the designing of Runge-Kutta methods, in particular, and the use 
of the Evolutionary Algorithm for any kind of optimization over an algebraic 
variety. To our knowledge this is the first time that algebraic geometry is used 
to state correctly that evolutionary algorithms have to be used to face a parti- 
cular optimization problem. Again we think this is the first time that algebraic 
geometry and evolutionary algorithms are used to tackle a numerical analysis 
problem. Further refinement of our evolutionary optimization algorithm will 
surely improve the solution of these important numerical analysis problem. 

A copy of the software can be obtained by sending an email to the authors. 
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02,1 
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62,2 


4 





11 
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4 


4 


9 


1 


2 


3 


5 


2 


9 


3 



Table 19: The Butcher tableau of the new solutions for a 3 levels Runge-Kutta of order 3; the error coefficients are of order 10 



02,1 


03,1 


03,2 


0.628 


0.076 


0.294 


0.736 


-0.082 


0.346 


0.706 


-0.198 


0.492 


0.755 


0.010 


0.233 


0.796 


-0.362 


0.566 


0.757 


-0.174 


0.416 


0.614 


0.076 


0.308 


0.695 


0.004 


0.299 


0.631 


-0.026 


0.395 



04,1 


04,2 


04,3 


0.921 


-0.330 


0.449 


0.834 


-0.237 


0.431 


0.952 


-0.322 


0.414 


0.753 


-0.185 


0.451 


0.914 


-0.279 


0.394 


0.857 


-0.246 


0.418 


0.947 


-0.352 


0.448 


0.846 


-0.259 


0.443 


0.977 


-0.362 


tO.434 



Wl 


W2 


W3 


0.159 


0.245 


0.433 


0.157 


0.286 


0.398 


0.149 


0.285 


0.412 


0.162 


0.283 


0.392 


0.145 


0.316 


0.390 


0.153 


0.298 


0.393 


0.159 


0.241 


0.438 


0.159 


0.269 


0.411 


0.155 


0.254 


0.432 



W4 


ei,o 


62.1 


62.2 


0.160 


6 


1 


4 


0.158 


8 


5 


4 


0.152 


1 


4 


4 


0.162 


6 


2 


4 


0.147 





3 


10 


0.154 


8 


15 


7 


0.160 


8 


6 


2 


0.160 


9 


1 


9 


0.157 


4 


15 


9 



63,0 


63,1 


63,2 


63,3 


11 


8 


4 


27 


13 


3 


9 


14 


1 


12 


20 


05 


2 


1 


5 


54 


3 





4 


10 


3 








5 


3 


9 


7 


14 


5 


1 


31 


4 


2 


5 


19 


7 



Table 20: The Butcher tableau of the new solutions for a 4 levels Runge-Kutta of order 4; the error coefficients are of order 10 
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